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We present numerical methods to solve the Israel-Stewart (IS) equations of causal relativistic 
dissipative fluid dynamics with bulk and shear viscosities. We then test these methods studying the 
Riemann problem in (1+1)- and (2+l)-dimensional geometry. The numerical schemes investigated 
here are applicable to realistic (3+l)-dimensional modeling of a relativistic dissipative fluid. 
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I. INTRODUCTION 

The interest in modeling the evolution of matter cre- 
ated in relativistic heavy-ion collisions with fluid dynam- 
ics has never ceased since the pioneering works by Lan- 
dau [l[ • Recent remarkable discoveries at the Relativistic 
Heavy Ion Collider (RHIC) at Brookhaven National Lab- 
oratory provide evidence for an almost "perfect" fluid- 
like behavior of the QCD matter created |2j. 

In a perfect, or ideal, fluid transport coefficients like 
bulk and shear viscosity and heat conductivity vanish. 
This is an idealized situation; in a real fluid one can show 
that there are lower bounds for these transport coeffi- 
cients, for instance using the uncertainty principle [3| or 
applying the AdS/CFT conjecture 0]. In order to decide 
how close to a perfect fluid the matter created at RHIC 
is, one has perform theoretical calculations in the frame- 
work of causal relativistic dissipative fluid dynamics. In 
this way, one may also be able to extract the numerical 
values for the bulk and shear viscosity coefficients from 
experimental measurements. 

Currently the most widely accepted and studied theory 
of relativistic dissipative fluid dynamics is due to Israel 
and Stewart [!, H, 0, H[ . This is the relativistic version of 
the pioneering work by Miiller [{J, [lOj. Although these 
theories have been developed in the 1970's, efforts to 
study and apply them to relativistic heavy-ion collisions 
have only started very recently [ll], [Tl] ■ This has been 
followed by an impressive number of studies in (1+1) 
dimensional 1X31 . Jl 4 . [l5l . [lal and (2+l)-dimensional ge- 
ometries 0, EI, El S IH H |M ll . 

In 3+1 dimensions, given arbitrary initial conditions 
and a general equation of state, the only way to solve 
the equations of relativistic fluid dynamics is by means 
of numerical methods. Any numerical method requires 
an algorithm that has to be tested in order to assess 
its validity for solving the underlying equations. Test- 
ing algorithms to solve relativistic dissipative fluid dy- 
namics is made difficult by the fact that there is only a 
rather limited number of test cases with analytical so- 
lutions. Reference [13] investigated sound propagation 
for the linearized IS equations. The algorithm of Ref. 
[ll[ was checked, for certain expansion scenarios, as to 
whether it correctly approaches the Navier-Stokes and 



ideal-fluid limits. So far, however, numerical algorithms 
to solve the IS equations have not been tested in situa- 
tions where shock discontinuities occur in the ideal-fluid 
limit. The present paper, in which we perform an ex- 
tensive study of the relativistic Riemann problem in 1+1 
and 2+1 dimensions, aims to fill this gap. 

In Sec. [U we provide a short review of IS theory of 
dissipative fluid dynamics. In Sec. Mil we formulate it 
in a form suitable for numerical implementation. This 
is followed in Sec. IIVI by an introductory presentation 
to the numerical methods which we use to solve the IS 
equations. In Sec. [V] we report results of solving the 
Riemann problem in 1+1 and 2+1 dimensions. Section 
IVII concludes this work with a summary of our results 
and an outlook. 



II. DISSIPATIVE FLUID DYNAMICS 
A. Units and definitions 

Throughout this work natural units, h = c — kg = 1, 
are used. Components of contravariant vectors and ten- 
sors in 4-dimensional space-time are denoted by upper 
indices, i.e., and A^ v . Greek indices take values 
from to 3 and Roman indices from 1 to 3. Covari- 
ant components, denoted by lower indices, are obtained 
by A v = g^vA^ , where g^ is the metric tensor, for which 
we use the (+,—,—,—) convention. If not stated other- 
wise the Einstein summation convention is used for both 
Greek and Roman indices. 

For an arbitrary contravariant four-vector A 11 the co- 
variant derivative is defined as 



At 



d a A» 



1 



(1) 



where 



\ g^ v {dpg a u + d a g uf3 - d v g a p) denotes the 
Christoffel symbol of the second kind and d a = d/dx a 
denotes the four-derivative. Similarly, the covariant 
derivative of covariant vectors is given by A M;Q = daA^ — 
T^ a A/3. For scalars the covariant derivative reduces to 
the ordinary four-derivative. The covariant derivative of 
second-rank contravariant tensors is 



d a A^ 



(2) 
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Vectors and tensors can be decomposed into parts par- 
allel and orthogonal to the four- velocity of matter u^, 
where u^u^ — 1. Using the transverse projection oper- 
ator A^ v = g^ v — u^v? where A^ v u v — 0, an arbitrary 
four- vector can be written as A 11 = u^UaA" + A^ a A a . 
The covariant derivative of an arbitrary tensor can be 
decomposed as 

A£-"" = u a DA^-^ + V a A^-^ , (3) 

where the convective time derivative D and the spatial 
gradient operator V Q are given by 

DA" 1 -"" = u^A^-^, (4) 
V a # 1 "'"" = A^A^- fln . (5) 

It is convenient to define the traceless and symmetric 
projection of a tensor field, which is orthogonal to u M . 
This is denoted by angular brackets (), 

A (H = Ia^A^(A q/3 + Af Ja ) - ±A^A af3 A al3 . (6) 

The covariant derivative of the four- velocity can be gen- 
erally decomposed as 

u U ;n = u^Duv + + -A^uO - w^v, (7) 

where the expansion rate 6, the shear tensor a^ Vl and the 
vorticity tensor tv^ are defined as 

e ee v^ = cy^ + r£X> (8) 

a** ee V<»u»> = l -A^A^{u a , p + u , a ) - °-A^ 
= \ - u»u a d a u" + d v u^ - u v u a d a u^) 
+ \ [a^vPy^ + A-Vr^) - 9 -A^, (9) 

wt = ^A MQ A^, {u a .f} — up. a ) (10) 
= - ( d v u» - d»u u + u»u a d a u u - u v u a d a u^ . 
where a^ v u u — and uu^ u u v = 0. 

B. The equations of causal relativistic dissipative 
fluid dynamics 

The basic quantities characterizing dissipative fluids 
are the net charge current and the e nerg y-momentum 
tensor T^ v . Following Refs. [H,ESE3,E1 these can be 
decomposed with respect to the fluid four-velocity u M as 

N" ee N? q + 5N» = nu" + V , (11) 
ee T£ + ST"" = eu^u v - [p + II) A"" 
+ W»u v + W u u" + tt^ , (12) 

where n = N^u^ is the net charge density and e = 
UfjT^Uv is the energy density in the local rest frame 



(LRF) , i.e., where = (1,0,0,0). The charge diffu- 
sion current is given by 5N^ ee V 1 = N U A^ V . The 
energy-momentum flow orthogonal to is given by 
ee A^Tupvfi . This quantity can be decomposed as 
W» = q^ + ( e +p)V^/n, where is the heat flow. The lo- 
cal isotropic pressure is denoted by p + II ee — ^A^T^", 
where p is the equilibrium pressure and II is the bulk 
viscous pressure measuring the deviation from the local 
equilibrium pressure. The shear stress tensor is defined 
as ee T^ v ' . This representation is completely gen- 
eral, valid in any coordinate system, and independent of 
the definition of the flow velocity. 

Usually, there are two typical choices used to define 
the flow velocity: either tied to the net charge flow when 
= (Eckart frame) or tied to the energy flow when 
W 1 * = (Landau frame) . We will use the latter definition 
in this work. 

Without conserved charges only Landau's definition of 
the flow velocity is appropriate. In this case the heat flow 
is q^ L = — (e +p)V fi /n. For net charge- free matter, q^ is 
not well defined, but also irrelevant for the discussion, so 
we set it to zero, q^ — — 0. 

When all dissipative quantities are zero, U M = W* = 
LT = it 1 - 1 " = 0, the decompositions (fTTj) and (fl2|) reduce to 
perfect fluid form, N 11 = N£ q = nu^ and T» v = T%» = 
eu^u v — p(e, n)A^ v . The LRF energy and charge densities 
are always fixed to their equilibrium values by the Landau 
matching conditions, i.e., n — n eq , and e = e eq . Then, 
the equilibrium pressure is given by the equation of state 
(EOS)p = p_(e,n) = -|A r T e f ; 

The equations of relativistic dissipative fluid dynamics 
follow from the covariant differentiation of the conserved 
charge four-current and the energy-momentum tensor, 

AT£ ee ±. 9li (JgN») = 0, (13) 

ee -L d, {y/gT^) + T^T^ = , (14) 

where g ee — det(<^„) is the negative determinant of the 
metric tensor. 

The non-equilibrium entropy current can be written as 

ee + 53* = su^ + , (15) 

where the entropy flux relative to u' 1 is = S V A^ V . 
The LRF entropy density is s = S^u^, where in general 
s < s eq (e,n). 

Following Refs. 0, H, [3, H, the phenomenological ex- 
tension of the entropy four-current by Israel and Stewart 
can be written without heat conductivity as 

S» EE SU» = S eq U» - (p Il 2 + (3 2 7T al3 7T a0 ) ^ , (16) 

where the coefficients /3o, /?2 are functions of e and n. 
Their exact value can be determined explicitly e.g. from 
kinetic theory. 

The requirement of non-decreasing entropy leads to re- 
laxation equations for the bulk pressure and shear stress 



3 



tensor. Here we also include the vorticity terms which 
follow from the kinetic-theory derivation, but we neglect 
the coupling between bulk and shear viscosity. Then, the 
IS equations @,[l6[ read 



DU 



— (Hjvs 
m 

1 



— (n 



n)- 



i l 



J 2 



J 3 



(17) 
(18) 



where rn = CPo denotes the relaxation time of the bulk 
viscous pressure and tv = 2r)[32 is the relaxation time 
of the shear stress tensor. The relativistic Navier-Stokes 
values are given by [HI, |26{ 



7T 



A r .S' 



-CO, 

2 W ^ . 



(19) 
(20) 



where £ > is the bulk viscosity coefficient and T) > 
is the shear viscosity coefficient. In Eqs. (fl7|) . (fl"8|) . we 
introduced the abbreviations 



Jo 






EE (tt^u" 


WW 

J 2 


"H 


WW 

J 3 


EE 27r^> 



ft 
T 



Dln^ 
T 



(21) 
(22) 
(23) 
(24) 



where we used tt^lv^ = 0. 

For the sake of simplicity, in our numerical studies pre- 
sented in the subsequent sections we assume a gas of 
massless Boltzmann particles without conserved charges. 
In this case, the equation of state is simply e = 3p 



and 



_ 3g 



T , where g is the number of degrees of 



freedom. The equilibrium entropy density is given by 
s eq = -^jT 3 . In this case, we can further simplify 
noting that the exact value of the thermodynamic inte- 
gral for massless Boltzmann gas is, ft = 3/(4p) [1, [30| . 
Therefore, it follows that £>/3 2 /ft = -De/e. Thus, 
Dln(ft/T) = -De/e - DT/T, where the temperature 
can be calculated from the EOS. The convective time 
derivative of the LRF energy density is given by energy 
conservation, De = — (e + P)utt — tt^u^, where the ef- 
fective pressure P is defined as P(e, n, II) = p(e, n) + II. 



C. General coordinate representation 

Here we give the relations between and N 11 in the 
calculational, or laboratory, frame and the LRF densities 
e, n, and the flow velocity v l . The natural frame of refer- 
ence is the laboratory frame. However, during the time 
evolution of the system, we have to extract the local ve- 
locity and the LRF densities from the laboratory frame 
quantities. These are needed because the EOS is given 
as a function of LRF densities, p = p(e, n). 



We can write the four-vector and tensor quantities 
given in Eqs. (fTTj) . (fT2"|) by specifying the four- velocity 
of the matter, — j(l,Vi) = j(l,v x ,v y ,v z ), where 



7 = (1 — v 2 ) 1 I 2 and v 
laboratory frame quantities take the form 

iV° 

T 00 = (r /' - -9 UU /' 

- 9 

= Vl T m + P{g m v l -g ^ ■ -"" 
T ij = (e + PrfvtVj-Pg 13 



n'yvi = ViN° , 

( e + P ) 7 2 „00, 
(e + P) 1 2 v l - q m P 



no 



V,.TT 
ij 



= Vl T°i+P(g 



. fl 'J)-« i7 r°3 + 



The 

(25) 
(26) 
(27) 

(28) 

(29) 



N° is the local charge density, N l is the local charge 
current in the direction i, i.e., the direction of the flow 
u L . The total energy density of the fluid is T 00 which in 
the LRF reduces to the (equilibrium) energy density e. 
By definition, T 0% denotes the energy flow in the direction 
of u\ while T l ° is the momentum density flux in the ith 
direction 1 . The remaining spatial part, T %3 , denotes the 
zth component of the momentum flowing in direction j. 

The LRF charge density and energy density are ob- 
tained from Eq. (J2SJ) and Eqs. l[2"7]). (J25J), respectively, 



n 
e 



N°( 

rpOO 



1 



2 )1 /2 j 

1 - Vi (T° l 



0i\ 



(30) 
(31) 



while Eq. |28J together with the above expressions leads 
to the expression for the velocity components, 



Pg 



0i 



yOO 



T 00 



Pg' 



oo ■ 



(32) 



In most cases of interest g 00 — 1 , and the metric of the 
space-time is diagonal. Therefore we can introduce a 
simplified notation which mimics the perfect fluid rela- 



y00 



T 00 



Mi 



T 



0i 



tions 29], R = nj, E 
where M = |M| = (M 2 + M 2 + Af 2 ) 1 / 2 . Thus, M is 
parallel to the velocity v, similarly as in the perfect fluid 
case. These quantities have to obey the physical con- 
straint M < E, in order to obtain meaningful solutions. 
Therefore, we can express the LRF charge density, en- 
ergy density, the absolute magnitude of the velocity, and 
the velocity components as 



n 

e 
v 

Vi 



R(l- vv) 1 / 2 . 
E- v-M, 
M/[E + P] , 
vMi/M . 



(33) 
(34) 
(35) 
(36) 



Substituting Eqs. ([33]). ([34]) into Eq. ((3SJ) we obtain the 
equation for the magnitude of the velocity, v. This can 



1 In standard units the flow of the energy density is cT 0t , while 
the flow of momentum density is c -1 T l °. 
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be solved by using a one-dimensional root search. There- 
after, use of Eq. (|36[) yields the individual velocity com- 
ponents and 7. Note that, in case of perfect fluids, this 
simplified treatment is practicable, however, in case of 
dissipative fluids, this may not always be possible. This 
is due to the fact that the vectors T ^ and ir 0)J " are not 
parallel to each other. Hence choosing other shear stress 
tensor components as independent variables, or in cases 
which take into account the heat flow, it is required to 
carry out a multidimensional root search to find the ve- 
locity [H, [13, [H, EH . 

For dissipative fluids the number of unknown variables 
increases by the introduction of the shear stress tensor 
and the bulk viscosity. The shear stress tensor is con- 
strained by the orthogonality condition tt^ v u v = 0, lead- 
ing to the following relations, 



n l0 U EE — TT^Uj 



nn o? ii l 

TT Uq — — TT Ui = IT J UjUi/Uo . 



(37) 
(38) 



One more independent relation follows from the trace of 
the shear stress tensor, tt^ 1 = 0, 



-TT Qi, 



(39) 



In general, using Eq. (|3T[) we can reduce the number of 
unknowns by three, and by two using Eqs. (|38|) and (|39|) . 
Thus we are left with five independent components of the 
shear stress tensor. However, for testing the numerical 
solutions it is preferable to calculate all shear stress ten- 
sor components directly using the relaxation equations, 
instead of using the orthogonality relations. We will re- 
turn to this matter later and provide examples. 



III. TEST PROBLEMS 

In this section we shall write the IS equations in various 
(1+1)- and (2+l)-dimensional geometries with Carte- 
sian or curvilinear coordinates. For the sake of com- 
pleteness, the (2+l)-dimensional boost-invariant and the 
(3+l)-dimensional IS equations in Cartesian as well as in 
(r, x, y, 77) coordinates are given in the Appendices. Here 
t is the longitudinal proper time and i] is the space-time 
rapidity. 



A. (1 + 1)— dimensional Cartesian coordinates 

In Cartesian coordinates, the metric tensor is g^ v 



.in' 



diag(l, — 1, — 1, — 1) and all Christoffel symbols 
vanish. The negative determinant of the metric is g = 1 . 
We assume that the system evolves along the z di- 
rection and that it is homogeneous in the transverse 
plane, such that the spatial derivatives in x and y di- 
rections vanish identically. The flow velocity of matter is 
u» = 7,(1, 0, 0, v z ) and 7, = (1 - v 2 z )- x ' 2 . 



N° 


= nj z 


} 






N z 


EE N°V 


Z 1 






1 


= (e + 


P)lz - 


p , ^00 _ 

r -\- 7T - — 


(e + P Z )7 Z - Pz 




ee (e + 




+ tt 0z = (e 


+ Pz)lzVz , 


rpXX 


EE TT XX 


+ p = - 


TT 

-2 +P < 




j^yy 


EE TT VV 


+ p = - 


TT 

-2 +P < 




rpZZ 


EE (e + 


P)llvl 


+ P + TT ZZ 


= (ev 2 z +P z h 2 z , 



The components of the energy-momentum tensor and 
charge current in the laboratory frame are 

(40) 
(41) 
(42) 
(43) 

(44) 

(45) 
(46) 

where the shear pressure, tt, is defined such that tt zz = 
7^7r. The orthogonality and tracelessness properties im- 
ply tt xx = ttVv = -tt/2, and tt 00 = v 2 j 2 tt, see Ref. [Hj]. 
From the orthogonality relation we obtain ir Qz = v z ir zz — 
v z j 2 tt. The effective pressure in the z direction is denoted 
by P z ee P + tt = p(e, n) +11 + tt. The remaining four- 
vector and tensor components vanish, N x — N y = and 

T 0x = T 0y = rpxy = Txz _ rpyz _ q This algQ meang 

that the corresponding shear stress tensor components 
vanish, tt 0x = tt^ = n xv = tt xz = n yz = 0. 

The LRF quantities and the velocity can be expressed 
in terms of the laboratory quantities, 

n = N°(l-v 2 z ) 1/2 , (47) 
e = T oa -v z T 0z , (48) 

Vz = roo + Pz ■ ( 49 ) 

The conservation equations follow from Eqs. (fT"3]) . (fT4|) , 

d t N° + d z (v z N°) = 0, (50) 
d t T 00 + d z (v z T 00 ) = -d z (v z P z ), (51) 
d t T Qz + d z (v z T 0z ) = -d z P z . (52) 

The relaxation equations for the bulk viscous and shear 
pressure follow from Eqs. (fT?)) . (p~5|) : 

7»3(n + w p = — (njvs - n) - 1 , (53) 

m 

-y z d t TT + j z v z d z TT = — (tt ns -Tx)-h, (54) 
T„ 

where If x — If x — 0. The Navier-Stokes values of the 
bulk viscous and shear pressure are 



-(Oz 

2-qa ■■ 



:V0z, 



(55) 
(56) 



V^u^ 1 denotes 



where 6 Z = d^vf 1 = da z + d z {^ z v z ) 
the expansion scalar and a = —2a xx = 8 z /3 is the shear 
stress. Furthermore Eqs. (HJ) and ([23]) with I 2 = -21% x , 
lead to 



In 



(57) 
(58) 
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B. (l+l)-dimensional cylindrical coordinates 

In the case of (l+l)-dimensional cylindrical coordi- 
nates, all quantities are functions of the time t and the 
radial coordinate r only. The flow velocity is given by 
= 7,-(l, u r , 0, 0), where 7 r = (1 — v 2 )^ 1 / 2 . The gradi- 
ent operator is 9 M = (<9 t ,<9 r ,0, 0). The terms containing 
drf, and d z vanish identically. 

The metric tensor transforms as g^ u — §§7r§fi7f7a,3, 
where x M = (t, r, <f>, z), — (t,x,y,z) and rjuu is 
the Cartesian metric. The spatial coordinates are r = 
,2 and 



y- 



arctan(j//x). The contravariant and 



covariant metric tensors are g^ v = diag(l, — 1, — 1/r , — 1) 



and g^ v = diag(l, -1, - 
ative determinant is g 
Christoffel symbols are r^ r = Tf^ = r _1 and 
The laboratory frame quantities are 



-1), respectively. The neg- 
2 . The only non- vanishing 



N° 
N r 

rpOO 

rprr 

r p<t>4> 

rpZZ 



(59) 
(60) 

P) 7 2 - P + 7T 00 = (e + P r ) 7 2 - P r , (61) 
P)^v r + n 0r = (e + P r ) 7 >r , (62) 



«7r ) 
jV°«, 

(e + 
(e + 

(e + P) 7 > 2 + P + 7r'' r 

(e + P r )7r«r + Pr , 

P . d>d> 
12 + ^ ' 

P + 7T ZZ , 



(63) 
(64) 
(65) 



where P r is the effective pressure in the radial direction 
defined below. All remaining vector and tensor compo- 
nents vanish identically, N* = N z = 0, = T 0z = 

T 4>r = T< f,z _ rprz _ q and ^00 = ^Oz = = rf>z = 

TT rz = 0. 

To reduce the number of unknowns we use the 
trans versality of the shear stress tensor, leading to ir 0r — 
v r ir rr , 7r 00 = v 2 Tr rr . The tracelessness condition gives 
7T 00 = 7r rr + r 2 ^ + n zz . The simplest solution is to 
choose tt^ and -k zz as the independent components of 
the shear stress tensor, since a Lorentz boost in radial 
direction does not affect these components. The remain- 
ing shear stress tensor components can be expressed by 
using these components, 



7T = 



-v rl 2 r {r 2 K^ + tt zz ), 
-v 2 rl 2 (r 2 7r^ + tt zz ). 



(66) 
(67) 
(68) 



The charge conservation equation and the equations 
of energy and momentum conservation follow from Eqs. 

da, m, 



d t N° + d r (v r N°) 
d t T 00 + d r (v r T 00 ) 



= —d r (v r P r ) 



- (v r T 
r 



00 



v r P r ) , 



(72) 
(73) 



d t T 0r + d r (v r T 0r ) = -d r P r 



2r 2 -K^ ~ tt zz ) . (74) 



Due to symmetry the right-hand side of Eq. (|74|) has 
to vanish at the origin. The relaxation equations follow 
from Eqs. m3), ([IS]), 

7r d t U + 7r v r d r U = i (Il N s - n) - la , (75) 
7 r d f 7r^ + lr v r d r TT rl " t ' 



— lV - 

Iff \ 

n lrVr 4,4, _ 



(76) 



lr d t ir zz + lr v r d r n zz = — (**> s -***)- I** ,(77) 

where the expansion scalar is 9 r — dt~f r + r~ 1 d r (r r ) r v r ) 
and if* = I{ z = 1^ = Pf = 0. Note that the convec- 
tive time derivative from Eq. (j4|) leads to an extra term 
for the ir** component, which was missed in Eq. (5.16) 
of Ref. 0. 

The shear stress tensor components are calculated from 
Eq. (O, hence the Navier-Stokes values for the bulk vis- 
cous pressure and shear stress tensor are, 



n 



NS 



'NS 



'NS 



= -ce r , 

= 2 Ti*++ - 2V ( 9r lrVr 
r 2 V 3 r 



2ry(T 2 



(78) 
(79) 

(80) 



Also note that the r~ 2 factor in 7rff s might cause prob- 
lems close to the origin. Hence it is preferable to 
rewrite the relaxation equation using the following vari- 
able: 7?^ = r 2 7r^. 

The term Iq and the relevant components of are 
given by 



The LRF charge density, energy density, and velocity are 
given as 



where P r = P 



n 


N°(l-v 2 r ) 1/2 , 


(69) 


i = 


T 00 - v r T 0r , 


(70) 


v r = 


rpQr 


(71) 


T 00 + P r ' 


TT rr 


P - rV^ - tt zz . 





To = 



T<P<P _ 



in 

2 



Pln^ 
T 



T 



0r + D\n 



(81) 
(82) 
(83) 



We also have the following relations between the cylin- 
drically symmetric and Cartesian systems (with similar 
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relations between the shear stress tensor components) 

T 0x = T Or cos0, (84) 



T°y = T 0r sin</>, 

T xx = T rr cos 2 (j> + r 2 T^ sin 2 , 

T yy = T rr sin 2 <j) + r 2 ^ cos 2 (j) , 

T xy = r T rr -r^T^) cos 6 sin 6. 



(85) 
(86) 
(87) 
(88) 

while T 00 and T zz remain unchanged. The inverse trans- 
formations are 

T or = T 0x cos(/> + T Vsin^, (89) 
T o<p = (T 0y cos <j>-T 0x sin (f)/r , (90) 
T rr = T xx cos 2 <l) + T xv sin(2(l>) + T yy sin 2 (j), (91) 
T** = [T xx sin 2 <j> - T xy sin(20) + T yy cos 2 </>] /r 2 , (92) 
T r<p = [(jyy -T xx ) sin (f> cos ct) + T xy cos(2(f>)}/r. (93) 

These relations will be used to compare the evolution of 
cylindrically symmetric and Cartesian systems. 

C. (2+1)— dimensional Cartesian coordinates 

For (2+l)-dimensional Cartesian coordinates, the co- 
variant derivative of any four-vector reduces to the stan- 
dard four-divergence, since all Christoffcl symbols van- 
ish. We assume that the system is homogeneous in 
the z direction and the velocity, as well as the deriva- 
tive in this direction, vanish. Hence the four-flow and 
four-gradient are function of (t, x, y) coordinates alone, 
thus u» = j±(l,v x ,v y ,0), <9 M = (d t ,d x ,dy,0), where 
7X = (1 - vl)- 1 / 2 and v x = (v 2 x + v 2 ,) 1 / 2 . 

The relevant laboratory frame quantities are 



N° 






(94) 


N x 




N°v x , 


(95) 


N y 




N Vy, 


(96) 


rpOO 




(e + Ph 2 ± - P + 7T 00 , 


(97) 


pOx 




(e + Py? ± v x + ir 0x , 








Vx T m + v x P-v x Tr ° + n 0x , 


(98) 


rpOy 




(e + Phlv y + ir 0y , 








VyT 00 +VyP-Vy7T 0a +n 0y , 


(99) 


rpXX 




(e + Phlv 2 x + P + 7T xx , 








v x T 0x + P- v x n 0x + ir xx , 


(100) 


pyy 




(e + P)ilv 2 y + P + n yy , 








v y T 0y + P- VyTr 0y + ir yy , 


(101) 


rpxy 




(e + P)jlv x v y + ir xy , 








v x T 0y - v x TT 0y + n xy , 








v y T 0x - v y ir 0x + ir xy , 


(102) 


rpZZ 




P + TT ZZ . 


(103) 



7r 0z — n xz — TT yz = 0. The LRF charge density and 
energy density are 

n = N° (1 - v 2 x ~ v 2 y ) 1/2 , (104) 
e = T 00 - TT 00 - v x (T 0x - tt 0x ) - Vy(T 0y - n 0y ) , (105) 

while the velocity components from Eq. (f32|) lead to 

rpOx ^OX 



x poo _ ^00 _|_ p 
pOy _ ^Oy 

u y = poo _ ^oo T p 



(106) 
(107) 



The velocity components can be calculated from the 
equations given above using a two-dimensional root 
search or via a one-dimensional root search using Eqs. 
(HH, ©I- 

Since we previously fixed and explicitly used some 
shear stress tensor components in the velocity calcula- 
tion, we choose to express the remaining components in 
terms of the former. The orthogonality relation ([57)1 and 
the tracelessness relation (|3"9"|) yield 



(108) 
(109) 
(110) 



Therefore, as function of the chosen independent vari- 
ables, 7r 00 , ir 0x , ir 0y , and ir zz , the other shear stress tensor 
components are 



7T 


= 1T XX V X 


+ TT XV V y , 


7T 


= TT Xy V x 


+ ir yy v y , 


7T 


= TT XX + 


„yy i „zz 



n xx = [v 2 y (n 00 - n zz ) + v x n 0x - v y ir 0y ] jv\ , (111) 
ttW = [^(tt 00 - ir zz ) - v x tt 0x - v y TT 0y ] jv\ , (112) 
n xy = [-VxVy(n 00 - n zz ) + v y 7T 0x + v x n 0y ] /v\ . (113) 



One may then check whether the remaining orthogonality 
relation, 



7T 00 = 7r° X Vx + !T 0y V, 



y ■■ 



(114) 



The remaining z-directed four-vector and tensor compo- 
nents vanish, i.e., N z = 0, T 0z = T xz — T yz — and 



is fulfilled. The above relations between the shear stress 
tensor components become unusable in the case that 
the velocity in the transverse direction approaches zero. 
Therefore, in our calculations we shall neglect the above 
simplifications and calculate all shear stress tensor com- 
ponents explicitly. 

Note that one can choose to select ir xx , iv yy , and ir xy 
as independent components, therefore 7r 00 , 7r 0:r , ir 0y , and 
7T ZZ are given by Eqs. (fl7)8j) - (fiT0)) and (fTT4|) , see Refs. 
[HI, [22j]. However, in this case the velocity iteration is 
two-dimensional which may become computationally as 
expensive as solving the additional transport equations. 

The conservation of net charge N°, energy T 00 , and 
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the momentum components T 0x and T 0y are 

d t N° + d x (v x N°) + d y (v y N°) = , (115) 
d t T 00 + d x (v x T 00 ) + d y (v y T 00 ) 



.00 i „0x\ 



(116) 



= -d x (v x P-V x W W +TT UX ) 
-d v (v y P-V y TT 00 +7T° V ), 

d t T 0x + d x (v x T 0x ) + d y (v y T 0x ) 
= -d x (P - v x n Qx + tt xx ) - d y (-v y TT 0x + n xy ) ,(117) 

d t T°y + d x (v x T 0y ) + d y (v y T°y) 

= -d x (-v x tt 0v + Tr xy ) - d y (P - v y TT 0y + n yy ) .(118) 

The relaxation equations for the bulk viscous pressure, 
II, and the components 7r 00 , -k 0x , Tr 0y , tt xx , n yv ,n zz , n xy of 
the shear stress tensor are 

7_i_<9 t n + yj_v x d x Tl + yxv y d y Tl 

= — (njvs - H) - Jo , (119) 

1± d t ir^ + 1± v x d x ir^ + 1± v y d y n^ 

= _1 _ tt^) _ l\» - 1$> - j£" , (120) 



where the Navier-Stokes values are IT 

T JVS 



NS 



-(9± and 



2r}a^ v , and the expansion scalar is 9± = <9t7_i_ 



d x ('j±v x ) + d y ("/j_v y ). 

The components of the shear tensor can be calculated 
from Eq. ([5]), which reduces to the following simple form 
a» v = i (d^u" - u^Du" + d^u^ - u v Du^) - in 
Cartesian coordinates. Hence, 

a 00 = d a± - 1± D 1± + ( 7 i - 1) ^ , (121) 



Jt.r 



r 0y _ 



2 [^(ti^x) - 9x7±] 

1 $ 

- [7±D(7x^) +7±^D7 ± ] +7i«x-g-, (122) 



[9t(7x%) - 9 y 7j_] 



2 [7x-D(7x« y ) + TJ-^-DT-i] + Tx^-jj" - ( 12 3) 



cr 2 ^ = -5 x (7j_t; x ) - 'yxVxDfrxVx) 
a yy = -d y (~f±v y ) - jj_v y D(~fj_v y ) 



(124) 
(125) 



The term If = (tt^u" + *k Xv u»)Du x leads to 

J? = 2 7± [it™D 1i _ - n 0x D( 7± v x ) 

- n 0y D( lx v y )] , (128) 
I°i x = 7± [(7r°°^+7r te )£>7x 

- (TT^ + Tr^Dfr^)] , (129) 
^ = 7± [(ir 00 v y + 7r 0y )D 1± 

- (TT° X V y +TT Xy )D( 1± V x ) 

- (tt°% + 7r ra )D(7 ± ^)] , (130) 
If' = 2 7± v x [ir 0x D 1± - n xx D( 1± v x ) 

- ir xy D( 1± v y )] , (131) 
If = 2 7 x% [7r^£>7_L - n xy D{ 7± v x ) 

- 7r yy D( 1± v y )} , 

I? = 7± [(^ + 7T% X )D7 ± (132) 

- {TT xy v y + n yy v x )D( 1± v y )] , (133) 

and Zf z = 0. The terms Iq and I^ v are again given by 
Eqs. (|21[) and (|23[) . respectively. Finally, the components 
of the term I^ v — k^u>\ + -k vX ^ x are explicitly given 
by 

If = 2 + 7r°^° y ) , (134) 

I 3 °* = Tr 0Q Lu x + Tr Qy Lu x y + Tr xx LU° x +iT xy u yl (135) 
I 3 ° y = 7r oa Lj y + 7r 0x Lj y x +7r xy Lj a x +n yy u; y , (136) 
If* = 2 (n 0x u x + Tr^w*) , (137) 
If = 2(n 0y u;y + n xy Lu y x ) , (138) 
I xy = tt 0x lo\ + n xx uj y x + n 0y uj x + n yy uj x (139) 



and /f = 



The vorticity tensor in 

Cartesian coordinates is given by, u)^ = 

| (<9j,m m — d^Uv + u^Du u — u^Du^), therefore the 
nonvanishing vorticity tensor components are 



u °x = \[9xl± + d t {j±v x )} 

+ g h± l 'x D l± - l±D(7±v x )] , 



(140) 



t x v 



[d X (j±Vy) + dy(^±V X )] 



2 h±v y D 1± - 7 ± D(7j.« v )] , (141) 



[~/±V x D(-~f±V y ) + 7xU v D(7±Uar)] 



^ = 2 - d x (-y±v y )] 



2 p J- 

3 ' 



where D = u^d^ = jj_d t + "f±v x d x + "f±v y d y . 



(126) + 2 ^v y D{"fj_v x ) - n/±v x D{ 7l _v y )] , (142) 

(127) where the vorticity tensor components satisfy the follow- 



ing relations, ui 
u> 0y = uiQy and uj : 



-uj 0x = uj 0x , = w y 



'„ = -u\ = -u xy = -u xy . 
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IV. NUMERICAL METHODS 

In this section we present in detail the numerical algo- 
rithm used to solve the equations of relativistic dissipa- 
tive fluid dynamics in (1+1)- and (2+l)-dimensional ge- 
ometries. In our case, this will be the SHArp and Smooth 
Transport Algorithm (SHASTA) [33j]. We also briefly 
discuss other schemes, and conclude with remarks on the 
numerical resolution and dissipative fluxes. 



A. One-dimensional implementation 

In (l+l)-dimensional systems the equations of charge 
and energy-momentum conservation, Eqs. |50|) . (f5l"j) . 
(|52[) . are of conservation type and can be generally writ- 
ten as 



d t U + d x (v x U) = S(t,x), 



(143) 



where U = U(t, x) is the conserved quantity, v x is the 
flow velocity in x direction, and S{t, x) is the source 
term. The relaxation equations (fB"3")) . are of con- 

vective type. These equations can be rearranged in con- 
servation form with an additional source term 13 li |32| , 



dtV-K + d x (v x Uir) = UAv x + S„(t, x) . 



(144) 



where U v is II or tt, and S v is the source term either from 

Eq. (SB]) or Eq. (54]) divided by 7 = (l - i^) _1/2 . 

To solve the above type of equations numerically, the 
original partial differential equations are replaced by ap- 
proximate algebraic difference equations and the values of 
U, v, and S are given at discrete grid points. The conser- 
vative, or primary, variable U(t, x) is replaced by its av- 
erage C/™ over the cell i at coordinate point X{, and at the 
discrete time step t n . The algorithms used in this work 
belong to the class of finite- volume methods where fluxes 
of the conserved quantity through the cell boundaries are 
calculated or approximated. This explicitly guarantees 
the conservation of the primary variable. The velocity 
and source terms are defined as a function of primary 
variables. Whenever source terms contain spatial deriva- 
tives, they are calculated by using second-order central 
differences, e.g. d x U? = {U? +1 - U^_ 1 )/(2Ax). Time 
derivatives in source terms are calculated using first-order 
backward differences, e.g. d t W? = (UJ 1 ' 1 - Uf)/At. 

Here we will give a brief presentation of our numeri- 
cal algorithm. Due to its simplicity, accuracy, and easy 
implementation for this study we choose the SHASTA 
[33| which was one of the first versions of Flux Corrected 
Transport (FCT) algorithms in the 1970's. Ever since, 
the FCT method has been extensively tested and refined 
for various studies, for example, the ETBFCT version by 
Boris 34], which also forms the basis for the LCPFCT 
algorith m 13511 . and the YDFCT algorithm by Toth and 
Odstrcil [361 ]. 

These explicit higher-order monotonic numerical meth- 
ods have been especially designed to work in the presence 



of strong gradients such as shocks. Typically low-order 
numerical schemes result in strong numerical diffusion 
due to the large truncation error, which tends to smooth 
out all the structures in the solution. Thus, low-order 
schemes are practically unusable unless unrealistically 
small grid sizes are used. Second-order schemes do not 
suffer from large numerical diffusion, but instead from a 
strong numerical dispersion, i.e., different Fourier modes 
propagate at different speeds. Especially in the presence 
of strong gradients like shock waves, numerical dispersion 
causes unphysical ripples in the solution, which eventu- 
ally invalidates the whole calculation. 

In the SHASTA this problem is solved by first calcu- 
lating a low-order solution which has a large numerical 
diffusion component. In the second step, as much diffu- 
sion as possible is removed from the low-order solution in 
such a way that no new maxima or minima are created, 
i.e., the monotonicity of the solution will be preserved. 
The remaining, or residual, diffusion of numerical origin 
is called numerical viscosity. In the FCT algorithms the 
numerical viscosity has both linear and non-linear contri- 
butions and therefore must be assessed separately for the 
problems at hand. This implicit numerical viscosity is, 
of course, different from the well-known explicit artificial 
viscosity of von Neumann, or Lax and Wendroff [37]. 

The low-order, or transported and diffused, solution in 
the explicit SHASTA method [33j is given by 



Ui 



1 



(Q 2 + A z - Q 2 _Ai 



(Q+-Q-)Ul l + AtS l 



Here, we defined 



Ai = U i+1 -Ui, 
_ _ l/(2A)Tt>? 



l/A±« ±1 -<) 



(145) 

(146) 
(147) 



where A = At/ Ax is the Courant number which in the 
SHASTA is restricted to values A < 1/2. The final time- 
advanced quantities are calculated by subtracting the so- 
called antidiffusion fluxes, A, from the transported and 
diffused solution such that 

U? +1 = U i -A i +A i -i, (148) 

where the flux-corrected antidiffusion flux is 



Ai = (j j max 



0, min ( CTjAi+i, \ Ai\,aiA 



i-l 



(149) 



Here, similarly as in Eq. (|146[) the difference of primary 
variables in adjacent cells is denoted by Aj = J7,+i — Ui, 
while the explicit antidiffusion flux is 



Ai = AadAJl 

cr, = sgn(Ai) . 



(150) 
(151) 



In the SHASTA, A ad = 1 is the default value of the 
so-called mask coefficient [38j]. This is a multiplicative 
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constant which can be set to lower values to reduce the 
amount of antidiffusion. 

Second-order accuracy in time is obtained by applying 
the SHASTA twice. First we calculate the velocity and 
source terms at time step n + 1/2. In the second step, 
these half-step velocity and source terms are used to cal- 
culate the final time- advanced quantity U" +1 . In a given 
cell, this can be summarized in formulas as 

yn+1/2 = V n ([fn w n At/2, Ax) , (152) 

U n+1 = U n (ll n , v n+1/2 , S n+1/2 , At, Ax) . (153) 

The relaxation equations are solved in a similar manner, 
however in this case the source terms actually depend on 
the primary variables, velocity field, and LRF quantities, 
therefore their values must be saved for full and half time 
steps. This requires much more memory compared to 
codes which solve ideal relativistic fluid dynamics. 



B. Multidimensional implementation 

The (2+l)-dimensional conservation equations are 
commonly written as 



d t U + d x (v x U) + dy{v y U) = S(t, x, y) 



(154) 



The cell-averaged conserved variable U(t, x, y) is denoted 
by Uf ■ • . A standard approach to solve such equations is 
to apply the dimensional or operator splitting method, 
which splits the original multidimensional equation into 
a sequence of (l+l)~dimensional problems [391 ]. 

A slightly different but more efficient approach is used 
in this work. We calculate the low-order transport solu- 
tion separately in the x and y directions by using the 
(l+l)-dimensional SHASTA (|145[) without the source 
term. Thus, the ^-transported quantity Ufj is given as 



Uf 



Ql 



1 r 
2 



(Q%y Af tJ ((>■ r a; , , 
(Q x + - ql) t/& , 

1/(2A*) T Ml, 



1/A* ± [(v x ) 



(155) 
(156) 



where Af d = U? +1}j - Ufa and \ x = Ax/ At < 0.5 is the 
Courant number in the x direction. A similar formula, 
with v x replaced by v y and all cell differences taken in y 
direction, holds for the y-transported quantity Uf- . The 
transported and diffused solution is then 



U, 



U x - 
^h3 



U v - 

w h3 



At S, 



h j 



(157) 



The advantage of this method is that it keeps the x — y 
symmetry of the system without the need to permute 
the directions in which the grid is updated. In this case 
it is also possible to implement a multidimensional flux 
correction in the FCT algorithm which avoids some nu- 
merical problems and leads to slightly smoother results 



compared to the dimensional splitting method for the 
same mask coefficient. To obtain second order accuracy, 
we use the method by DeVore 14011 , which is an improved 
version of Zalesak's method [4JJ. The full solution is 
given by 



U 



71+1 



AX 
A %3 



A 



(158) 



where the A's are the limited antidiffusion fluxes given 
in Eqs. (fTTTj) and (fT72| below. 

As in the (l+l)-dimensional case the antidiffusion 
fluxes in x and y directions are given by 



.4 



1,3 



(159) 
(160) 



where A x d ,A^ d are the antidiffusive mask coefficients, 
similarly to the (l+l)-dimensional case. Furthermore, 



A' 



U, 



i+1,3 



■u, 



1,3- 



(161) 
(162) 



In the DeVore scheme, the antidiffusion fluxes in x and 
y directions are first limited as in the (l+l)-dimensional 
case, 



A x 

A i,3 



a X j max [0, 

mm (<,A* +1J , lAfjlafjAU,^} , (163) 
<j\ j max [0, 

min (< J .A|' J . +1 , Al^^ :j ,)" , (164) 



where af j = sgn(Af J ) and afj = sgn(A^). Note that 
this additional step was introduced by DeVore into the 



multidimensional flux limiting algorithm by Zalesak. 

The allowed values for Uf't 1 after the antidiffusion 
stage are between 



jjmin 
U i,3 



min (Oij-i, Ui-ij, U i:j , U i+ ij, t/jj+ij , (165) 

(Uij-i,Ui-ij,Uij,Ui+ij,Uij+ij . (166) 

The total incoming and outgoing antidiffusive fluxes in 
cell are calculated as 

A™j = max ^0, A^_ ld ^ ~ min (o, Af A 

+ max (0, A^-.J - min (o, , (167) 

A?f = max (0, A X J\ - min (o, A?_ hj 

+ max (0, A\ \ - min (o, • (168) 

This information is then used to determine the fractions 
of the incoming and outgoing fluxes, 



pin _ 

i,3 



rpout 
i,3 



U i,j — U i,3 



(Uij - ' 77") 



/A 



1,3 



Mout 



(169) 
(170) 
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which is subsequently limited so that it creates no un- 
dershoot or overshoot in the cell it is leaving or entering. 
Thus, the new antidiffusive fluxes are given as 



-4? 



and 



.mm(l,F™ hj ,F?f), if Af 3 >0, 

« i min{ !./;</;. /;■:<;. ; 



3 - (171) 
if M 4 < 0, y ' 



p 7v fmin(l,i^ +1 ,i^*), if 4,>0, 

A iJ X j min ( 1; Fin^ pout_J i ^ < Q _ U ^ 

This (2+l)-dimensional numerical scheme can be gener- 
alized to (3+1) dimensions by extending the method to 
another spatial direction and repeating the above steps. 



C. Other numerical schemes 

Computational fluid dynamics (CDF) is a constantly 
growing held of research. There is a vast amount of 
methods which have been designed to solve the special 
relativistic fluid-dynamical equations in the perfect fluid 
limit, see Ref. [42j and references therein. 

In applications to relativistic heavy-ion collision the 
FCT-SHASTA and RHLLE methods have been system- 
atically explored for various test problems and shown 
to give excellent agreement [H, |3 EH . This is one of 
the reasons why we have chosen the SHASTA for our 
study. There are other well-known methods widely used 
in astrophysics and heavy-ion physics, such as Smoothcd- 
Particle Hydrodynamics (SPH) [46[ which has been re- 
cently extended to dissipative fluids [47[ , or the Particle- 
In-Cell (PIC) method [48| , but since they are completely 
different from finite-volume schemes we will not go into 
details. 

Other methods of interest such as High-Resolution 
Shock-Capturing (HRSC) methods based on the exact 
or approximate Riemann solution proved to be superior 
to SHASTA [H, |H, Hi| . However, in case of dissipative 
fluids such methods become difficult to apply due to the 
fact that there are no known analytic or approximate so- 
lutions for the Riemann problem. Recently, new methods 
have been developed to solve the hyperbolic equations of 
conservation or relaxation type which sidestep the need 
of Riemann solvers but have an accuracy comparable to 
HRSC schemes. These are new High-Resolution Central 
Schemes (HRCS) improving on the Lax-Friedrichs cen- 
tral scheme [H(J. The most important of these are the 
Nessayu-Tadmor (NT) and the Kurganov-Tadmor 
(KT) 52| schemes, see Ref. [53| for a collection of refer- 
ences. 

The KT scheme improves upon the NT scheme using 
information about the local propagation of speeds, which 
becomes problematic to evaluate for the IS equations. 
However, it gives excellent results for perfect fluids [54| . 
An important extension of the NT scheme was made by 
Pareschi [55| to describe both the stiff and unstiff regions 



of hyperbolic relaxation equations such as the IS equa- 
tions or equations of Ottinger-Grmela type (56|. In the 
latter case, this method has been shown to provide ro- 
bust results and excellent agreement between the (1+1) 
and (2+l)-dimensional cases [57j • Following this work 
we also made use of both the NT and KT schemes and 
compared them with SHASTA for the (l+l)-dimensional 
evolution of a perfect fluid. The results are very robust 
and agree very well. Therefore, without much more effi- 
cient methods at hand we simply choose to solve the IS 
equations with the SHASTA. 



D. Remarks on numerical resolution 

Fluid dynamics is a theory which is valid on time and 
length scales which are larger than the underlying micro- 
scopic time and length scales. In solving the equations of 
fluid dynamics numerically, we should be able to resolve 
all relevant time and length scales in the problem. In 
practice this means that the grid spacing Ax and time 
step At should be smaller than any of these scales. In 
perfect fluid dynamics, or in the Navier-Stokes theory, 
all scales are macroscopic, i.e., they are inversely propor- 
tional to the gradients of the fluid-dynamical variables 
like flow field and densities. Thus it is sufficient to have 
a numerical resolution that correctly resolves the macro- 
scopic structures. 

In the IS theory we also need to solve the relaxation 
equations for the dissipative currents. In this case the 
relevant time scale to be resolved is the relaxation time 
tr, which is of the order of the mean time between the 
collision of particles. Thus, the time step should be cho- 
sen such that At <C tr. If tr is much smaller than the 
macroscopic scales, this might require very high resolu- 
tion and therefore lead to very demanding calculations. 
However, in modeling heavy ion collisions, an application 
which we mainly have in mind, scale separation by sev- 
eral orders of magnitude is not expected throughout the 
whole fluid-dynamical evolution. 

There exist specialized methods 55] to solve the equa- 
tions in stiff regions. However, we do not consider these 
methods here, but simply choose sufficiently high res- 
olution to resolve both the macroscopic and relaxation 
time scales. Therefore, we solve simultaneously both the 
conservation and the relaxation equations with the same 
numerical resolution and scheme. 



E. Remarks on dissipative fluxes 

In relativistic dissipative fluid dynamics, the compo- 
nents of T A " / and cannot take arbitrary values. Obvi- 
ous physical constraints are that the LRF energy density 
must be positive semi-definite and the velocity must be 
bounded from above by the speed of light, i.e., e > and 
v < 1. Another constraint follows from the equation for 
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energy conservation, 

d^T " = d t T 00 + V • (vT 00 ) = , (173) 

where v % = T (H /T 00 . In order to have causal propagation 
of energy, we have to require that |v| < 1, i.e., 

V-T 0l T 0l < T 00 . (174) 

For perfect fluids, because of Eq. (|49|) . this condition 
guarantees both e > and v < 1, provided that the pres- 
sure is positive. However, in dissipative fluid dynamics 
this is not necessarily true, since the condition (|174p is 
sufficient only if the effective pressure is positive. For 
example, neglecting the shear pressure this leads to the 
condition II > — p for the bulk viscous pressure. 

The IS theory does not itself restrict the values of the 
dissipative quantities. In principle any value of shear 
and bulk pressure can be used, e.g. as an initial condi- 
tion. However, the applicability of the theory requires 
that the dissipative currents give sufficiently small cor- 
rections to the equilibrium quantities. For the shear and 
bulk pressure this requirement can be stated as 

|II| < Cp, (175) 
KH < C\Tg\, (176) 

where C is a constant of order, but smaller than, one. If 
these conditions are not satisfied, fluid dynamics is not 
expected to give a reasonable description of the space- 
time evolution of the system and the numerical calcu- 
lation can become unstable 32]. To protect the code 
from these numerical instabilities the conditions (|175p 
and (|176p are always enforced. This means that after 
each time step the above conditions are checked and 
tt^ v and n are adjusted accordingly. We note that the 
above conditions may be enforced before the velocity root 
search, in which case we have to compare to the values of 
T£q andp from the previous time step. Alternatively, one 
can apply these limiters inside the root search algorithm. 
In this case the limiters are applied simultaneously with 
solving for the LRF densities and the velocity. In every 
iteration of the velocity root search the shear and bulk 
viscous pressure are compared to the values of and p 
at the current time level. This guarantees that the con- 
ditions are always fulfilled, but the drawback is that this 
is computationally more expensive. In situations where 
we expect fluid dynamics to give a reasonable descrip- 
tion these conditions need not be enforced. However, if 
they are violated only in small regions of space-time, i.e., 
few cells or few time steps, enforcing the inequalities can 
prevent these regions to invalidate the whole calculation. 
Naturally, if the inequalities are violated in large regions 
of space-time, it signals that fluid dynamics is no longer 
a valid theory for such situations. 

V. RESULTS OF COMPARISONS 

In this section we apply the different numerical 
schemes described above to the relativistic Riemann 



problem in (1+1) and (2+1) dimensions. In (1+1) di- 
mensions the Riemann problem is analytically solvable 
for perfect fluids. Thus, it provides an important test 
case to compare the performance and accuracy of differ- 
ent numerical algorithms. 

Unfortunately, analytic solutions for the one- 
dimensional viscous Riemann problem are, to the best 
of our knowledge, not known. However, this type of 
one-dimensional test was performed previously: our 
fluid-dynamical calculations with non-zero viscosity 
were shown to give good agreement with kinetic theory 
simulations using the Boltzmann Approach to Multi- 
Parton Scatterings (BAMPS) [58[ parton cascade code 
[59L |60| . The purpose of our tests here are to show 
that a more complex (2+l)-dimensional code can, with 
similar initial conditions, remarkably well reproduce our 
earlier results for (1+1) dimensions. This confirms that 
the numerical method produces correct answers in these 
test scenarios, and gives us confidence that it can be 
successfully used to study phenomena where dissipation 
plays an important role. 

We shall proceed as follows: First, the Riemann prob- 
lem is briefly introduced and its analytic solution in 
(1+1) dimensions is compared with numerical solutions 
in the perfect fluid limit. Here we compare the SHASTA, 
the NT, and the KT numerical schemes. They all give 
comparable results and can reproduce the analytic re- 
sults with sufficiently good numerical resolution. This 
gives confidence that any of the schemes forms a good 
basis to extend the calculation to multidimensional prob- 
lems as well as to non-zero viscosity. In this work these 
extensions are made by using the SHASTA. 

Second, the numerical solutions for the (1+1)- 
dimensional Riemann problem with non-zero shear vis- 
cosity are shown. We compare results from the 
(2+l)-dimensional code to the results from the (1+1)- 
dimensional code and show that both codes yield, to good 
accuracy, the same results. 

Finally, the numerical solutions of the (2+1)- 
dimensional, azimuthally symmetric Riemann problem 
with non-zero shear viscosity are studied. We compare 
the results from the (l+l)-dimcnsional azimuthally sym- 
metric code to the results from the (2+l)-dimensional 
code. Again, these calculations are in excellent agree- 
ment with each other. 



A. The Riemann problem 

The initial setup for the (l+l)-dimensional Riemann 
problem consists of two states with constant pressure, po 
and pi, separated by a membrane at z — 0. The matter 
is initially at rest on both sides and homogeneous in the 
transverse directions. After the membrane is removed, in 
thermodynamically normal matter [44| there is a shock 
wave traveling into the region with lower pressure, and a 
rarefaction fan into the region with larger pressure. The 
interface between the two regions moves at a constant 
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velocity and is called the shock plateau. In dissipative 
fluids due to non-zero viscosity the initial sharp disconti- 
nuity will be smeared out and the quantities will change 
smoothly rather than discontinuously. 

In numerical calculations, unless stated otherwise, we 
have fixed the parameters as follows: The local Courant 
number is X x = X y = 0.4, and the comparison is made 
at t = 4 fm. The cell sizes Ax, Ay, and the antidiffusion 
coefficients, A x ad , A v ad are specified separately in all cases. 

The energy density in local equilibrium is given by 
e = ff T 4 where g = 16 is the number of degrees of 
freedom. Therefore, on the left and right-hand side of 
the initial discontinuity the energy densities correspond 
to the following temperatures: on the left T = 0.4 GeV 
and on the right T4 = 0.2 GeV. The bulk viscosity to 
entropy density ratio £/s and the shear viscosity to en- 
tropy density ratio 77/ s are taken as a constant, where 
the entropy density s = s eq is fixed to its equilibrium 
value, s eq = ||T 3 . In all test cases we start from local 
thermal equilibrium, i.e., initially — 0. We show only 
results with shear viscosity, but we have tested that we 
get similar results with non-zero bulk viscous pressure. 



B. Comparing different methods in perfect fluid 
dynamics 

The first test compares how well the different numer- 
ical methods can reproduce the analytic Riemann solu- 
tion [43[ in the perfect-fluid limit. The left panel of Fig. 
[T] shows the velocity v, the LRF energy density e, and the 
expansion rate 9 calculated with the SHASTA, the NT, 
and the KT schemes compared with the analytic solution. 
The numerical calculations in the figure are made with 
cell size Ax = 0.1 fm and At = 0.04 fm/c. We used the 
non-staggered version of the HRCS schemes with a min- 
mod limiter (6 — 2) which ensures that no local extrema 
are introduced, see Eq. (4.9) in Ref. [52]. 

All algorithms reproduce the analytic solution with 
nearly the same accuracy and numerical artefacts. In 
particular, all methods show long-wavelength oscillations 
which are best visible in the expansion rate, in the region 
between the rarefaction tail and the shock wave. The 
HRCS calculations show a somewhat larger overshoot 
for the velocity at the contact discontinuity as well as 
a more diffused shock front compared to SHASTA with 
A ad = 1.0. 

We also compared the above SHASTA result with a 
calculation with a reduced mask coefficient A a( i = 0.8, 
shown in the right panel of Fig. [1] This reduction 
strongly suppressed the unphysical oscillations in the nu- 
merical solution, but leads also to more diffusive pro- 
files. Furthermore, with the standard mask coefficient 
we have used the viscous SHASTA (vSHASTA) 2 solver 



with a small physical viscosity 77/s = 0.01 and A a d = 1.0. 
This very closely reproduces the A a d — 0.8 results with 
77/s = 0, especially at the smooth parts of the solution. 
Therefore, albeit small discrepancies exist at the shock 
font, we can conclude that our numerical solutions with 
the reduced antidiffusion mask coefficient have an addi- 
tional numerical viscosity corresponding to 77/s ss 0.01 
compared to the A a d = 1.0 case. 

Since all numerical calculations only approximate the 
exact solution, there is always some residual numerical 
viscosity in the solution. In fact, some amount of numer- 
ical viscosity is required to stabilize the solution. How- 
ever, this residual numerical viscosity can be made arbi- 
trary small by increasing the resolution. This is demon- 
strated in Figs.[2Ja), (c), (e), where all numerical algo- 
rithms considered reproduce the analytic solution almost 
perfectly with a cell size of Ax = 0.01 fm and At — 0.004 
fm/c. Also, the additional numerical viscosity result- 
ing from the reduction of the mask coefficient A a d scales 
approximately with the cell size for a constant Courant 
number. This is demonstrated in the right panel of Fig. 
[3J where we found that the additional numerical viscos- 
ity corresponds to 77/s » 0.001. We have checked that 
we get similar results also with other initial temperature 
ratios. 



C. Comparison between the one- and 
two-dimensional solutions 

The next numerical tests consist of comparing the 
(1+1) -dimensional solution to the (2+l)-dimensional so- 
lution of the one-dimensional Riemann problem in Carte- 
sian coordinates. The one-dimensional Riemann problem 
can be initialized on a two-dimensional grid in several 
different ways. We study here two different initializa- 
tions. In the first case, the initial discontinuity is along 
the y axis, i.e., on the x = plane. In the second case 
we place the discontinuity on the y = — x plane. These 
two cases are compared to the (l+l)-dimensional cal- 
culation. Here both one- and two-dimensional calcu- 
lations are done using the vSHASTA algorithm, with 
A a d = 0.8, grid size Ax — 0.2 fm and non-zero shear 
viscosity 77/s = 0.1 in all cases. 

In the simple one-dimensional formulation, there are 
only two dissipative quantities to propagate, Hi and 
7i"i , while the other shear stress tensor components are 
straightforward to calculate. In the two-dimensional 
setup we always propagate all non-vanishing dissipative 
tensor components, n 2 , ir® , n® 1 , , 7rf x , flf 9 ' 1 ^2* ■ 
Because there is only one independent shear stress com- 
ponent 7r in the one-dimensional Riemann problem, any 
of the non- vanishing shear stress components in the two- 



2 Our abbreviation only specifies that next to the conservation 



equations we also solve the relaxation equations of the physical 
viscosity using SHASTA. 
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FIG. 1: (Color online) The analytic (thin line) and numerical solutions of the relativistic Riemann problem on a grid with 
N x — 100 cells with Ax — 0.1 fm, after Nt = 100 time steps at t = 4 fm/c. (a) the collective flow velocity of matter, v, 
calculated with the SHASTA (continuous line), and the NT (dashed line) and KT (dotted line) algorithms, (b) The velocity, 
v, calculated with SHASTA using a mask coefficient A a d = 1.0 (continuous line), A a d = 0.8 (dashed line), and vSHASTA with 
A a d = 1.0 and rj/s = 0.01 (dotted line). Similarly, the LRF energy density, e, and the invariant expansion rate, 9, are shown 
in the panels (c), (d), and (e), (f), respectively. 



dimensional calculations can be used to extract n. The 
simplest possibility is to use tt = —2tt zz , because tt zz is 
independent of the orientation of the initial state in the 
(x, y)-plane. 

The result of the comparison between the one- and the 



two-dimensional calculations is shown in Fig. [3l where 
we compare the velocity v, the LRF energy density e, the 
expansion rate 9, and the shear pressure w. The velocity 
v = v z in the one-dimensional calculation, while v = v x 

or v — t/v'i + v?. in the two-dimensional cases. In the 
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FIG. 2: (Color online) As in Fig. [fl except for N x = 1000 cells with Ax = 0.01 fm, after Nt = 1000 time steps at t = 4 
fm/c. Analogously, the shear viscosity to entropy density ratio in the vSHASTA calculations shown in panels (b), (d), (f) is 
r)/s = 0.001. 



two-dimensional calculations the quantities are plotted 
along the x axis when the initial discontinuity is at x — 0, 
or along the y — x line when the discontinuity is in the 
y = —x plane. 

When the initial discontinuity is in the x = plane the 
two-dimensional SHASTA reduces essentially to the one- 
dimensional one. This is because there are no gradients 
in the y direction and v y = 0. Therefore, in this case 



we expect very good agreement between the one- and 
two-dimensional calculations. This is confirmed in Fig. 
[31 where the two-dimensional calculation (dashed line) is 
basically on top of the one-dimensional calculation (solid 
line) . 

When the initial discontinuity is along the y = —x 
plane, there are gradients in both x and y directions 
and both velocity components v x and v y are non-zero. 
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FIG. 3: (Color online) The numerical solution of the relativistic Riemann problem with rj/s = 0.1 on a symmetric grid with 
N x — N y = 200 cells with Ax = 0.2 fm, after N t = 50 time-steps at t = 4 fm/c. In all figures, the full line shows the 
one-dimensional evolution of matter. The dashed line shows the two-dimensional solution in x direction, while the dotted line 
shows the solution in the diagonal x — y direction, (a) The collective flow velocity of matter, v, (b) the LRF energy density, 
e, (c) the invariant expansion rate 9, and (d) the shear viscous pressure, n. 



This calculation is shown as dotted line in Fig. [3] The 
agreement with the (1+1) -dimensional results is still 
very good, although the two-dimensional algorithm gives 
somewhat sharper profiles in the shock region. 

The next test compares the (1+1) -dimensional solu- 
tion in cylindrically symmetric coordinates from Sec. 
IIIIBI against the two-dimensional solution in Cartesian 
coordinates with cylindrically symmetric initial condi- 
tions. This tests how well the two-dimensional system 
keeps its symmetry in time and performs compared to 
the one-dimensional counterpart. 

The initial discontinuity lies on a circle with radius 
ro = 5 fm, with a cell size of Ar = 0.2 fm in both 
cases. The velocity and position in the one-dimensional 
case is v — v r and x = r, while in the two-dimensional 

case v = ^Jv 2 , + v 2 and r = \J x 2 + y 2 . The first two- 
dimensional result compares the evolution of the system 
along the x axis, i.e., y = 0, while the second one does 
this along the diagonal, x = y. These are plotted with 
dashed and dotted lines, respectively, against the one- 



dimensional solution (solid line) in Fig. |U The other 
plots show the expansion rate, and the shear stress ten- 
sor components, ir zz , r 2 ir^, and ir 0r as calculated from 
the different equations in Sees. IIIIBl and UlI CI 

Similarly as before, the results are nearly the same, 
however, differences in the diagonal direction are visible 
and more pronounced than along the coordinate axis, due 
to the finite resolution. The agreement will obviously get 
better by decreasing the cell size and time step. 



VI. CONCLUSIONS 

In this paper, we have studied numerical algorithms 
to solve the IS theory for relativistic dissipative fluid 
dynamics. First, we briefly reviewed the IS theory 
and wrote the IS equations for (1+1)- and (2+1)- 
dimensional systems in Cartesian coordinates, and for 
(l+l)-dimensional azimuthally symmetric systems in 
cylindrical coordinates. For the sake of completeness we 
present the (3+l)-dimensional equations in Cartesian co- 
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FIG. 4: (Color online) The numerical solution of the relativistic Riemann problem in cylindrical geometry, with rj/s — 0.1 on 
a symmetric grid with N x = N y — 200 cells with Ax = 0.2 fm, after Nt ~ 50 time-steps at t = 4 fm/c. In all figures, the full 
line shows the one-dimensional evolution of matter. The dashed line shows the two-dimensional solution in the x direction, 
while the dotted line shows the solution in the diagonal x = y direction, (a) The collective flow velocity of matter, v, (b) the 



LRF energy density, e, and (c) the invariant expansion rate, 0. The shear stress components n z 
panels (d), (e), and (f), respectively. 



and 7r are shown in 



ordinates and the (2+l)-dimensional boost-invariant and 
(3+l)-dimensional equations in (r, x, y, rf) coordinates in 
the Appendices. We also gave a detailed introduction 
to the FCT-SHASTA method for one- and multidimen- 
sional applications, together with a brief discussion on 



the HRCS methods NT and KT. We also discussed re- 
lationship between microscopic and macroscopic scales, 
as well as physical limitations for the components of the 
energy-momentum tensor. 

In our first numerical comparison we solved the (1+1)- 
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dimensional Riemann problem in the perfect-fluid limit. 
This problem has an analytic solution which allowed us 
to make a definite comparison of performance and ac- 
curacy of the different numerical algorithms. For this 
problem all the algorithms considered here, i.e., the NT, 
KT, and SHASTA methods, gave very similar results. 
All of them could reproduce the analytic solution with a 
very high precision with sufficiently high numerical reso- 
lution. Moreover, with the same resolution the accuracy 
of the methods was found to be similar, i.e., none of 
them showed significantly faster convergence to the ana- 
lytic solution when the grid spacing was decreased. For 
this reason we have chosen SHASTA for all the other 
geometries as well as for all calculations with non-zero 
viscosity. 

We further studied the effect of the mask coefficient 
A a d in the SHASTA. This numerical parameter controls 
the amount of numerical diffusion in the algorithm. It 
was found that a reduction of the coefficient by 20% from 
the default value smoothens unphysical sharp structures 
in the solution, especially in the expansion rate, and at 
the same time only increases the numerical viscosity by 
a small amount. 

In the case of non-zero viscosity, the analytic solution 
to the Riemann problem is not known. However, we 
have demonstrated earlier that our (l+l)-dimensional 
code is in good agreement with kinetic-theory calcula- 
tions [m [60]. In this work we have first applied both 
(1+1)- and (2+l)-dimensional Cartesian implementa- 
tions of the code to the same (l+l)-dimensional Rie- 
mann problem. In this case we have chosen a non-zero 
shear viscosity, rj/s = 0.1. If the discontinuity in the 
initial energy-density profile was chosen to be along one 
of the coordinate axes, perfect agreement between the 
one- and the two-dimensional codes was found. In the 
case where the initial discontinuity was chosen to be 
along the y — —x plane, a slight difference between the 
two codes near the shock front was found. The results 
were shown for a rather large grid spacing Ax = 0.2 
fm; the agreement was found to improve significantly 
for smaller grid spacing. A similar comparison between 
the (l+l)-dimensional solution in cylindrical coordinates 
versus the (2+l)-dimensional solution in Cartesian co- 
ordinates with cylindrically symmetric initial condition 
confirmed that our method works well also for problems 
in more than one spatial dimension. 

In this work, we have demonstrated the applicability 
of FCT-SHASTA to solve the conservation equations of 
causal relativistic dissipative fluid dynamics simultane- 
ously with relaxation transport equations. In the fu- 
ture, we intend to extend this method to full (3+1)- 
dimensional geometries. We plan a detailed comparison 
with calculations done in the framework of kinetic the- 
ory [6l| , as well as studies of collective flow in relativistic 
heavy-ion collisions. 
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APPENDIX A: (2+l)-DIMENSIONAL 
BOOST-INVARIANT EXPANSION 

Because it is very important for modeling ultrarel- 
ativistic heavy-ion collisions, we discuss the (2+1)- 
dimensional boost-invariant equations of motion. The 
metric tensors are g**" = diag(l, — 1, — 1, — 1/r 2 ) and 
9iiu = diag(l, — 1, — 1, — r 2 ), leading to g = t 2 , where 
t = (t 2 — z 2 ) -1 / 2 is the longitudinal proper time and 
i] = l/21n [(t + z)/(t — z)] is the space-time rapidity 
(which is not to be confused with the shear viscosity co- 
efficient). The only nonvanishing Christoffcl symbols are 
P> = T'' = t -1 and T T = t 

The equations of relativistic dissipative fluid dynamics 
can be easily derived from the results in Cartesian coordi- 
nates, cf. Sec. IIII CI In order to obtain the equations for 
the boost-invariant case, the indices for time t and spatial 
z direction have to be replaced by r and rj in the four- 
vector and tensor components, (0,2:, y, z) — > (T, x,y,r)). 
Therefore, we easily find that all laboratory frame quan- 
tities can be written in the same way as in Eqs. (IM)) - 
(I102p , with the exception of the term in Eq. (|103|) , which 
becomes, T zz -> T m = P/t 2 + tt''''. This means that 
N n = 0, T T, > = T xr > = P'' = 0, and tt™> = ir xr > = = 
0. 

The LRF charge density, energy density, and velocity 
are calculated the same way as in Sec. IIII CI The charge 
conservation equation is 

d T N T + d x (v x N T ) + d y (v y N T ) = ——N T . (Al) 

T 

The energy conservation equation follows from 
d T T TT + d x (v x T TT ) + d v (v y T TT ) 

= -d x (v x P - V X TT TT + 7T TX ) - dy (VyP - VyTT TT + lT T y) 

_I (T TT + P + tV) . (A2) 
The momentum conservation equations follow from 
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d T T TX + d x (v x T TX ) + d y (v y T TX ) = --T™ 
- - 8 X (P- v x tt tx + n xx )- 8 V (-v y n TX +tt xv ), (A3) 



d T T Ty + d x (v x T Ty ) + d y (v y T Ty 



-_rpry 



- d x (-v x ir TV + ir xy ) -d y {P- v y Tr Ty + ir yy ) . (A4) 

The use of boost-invariant coordinates affects the expan- 
sion rate, ± = 7j_/t + d T ^± + d x (j±v x ) + dy(j±v y ), 
and the a zz component of the shear tensor, which is 



replaced by a 



vv = T-- 2 



9j_ /3 — jj_ /t) . The convective 



time derivative D from Sec. IIII CI becomes D = j±d T 
j±v x d x + r y±v y dy. The relaxation equations are the same 
as in Cartesian coordinates except for the replacement 
ir zz — ► tt vv which due to a nonvanishing Christoffel sym- 
bol includes a new term, 27t ?, ' ? 7_l/t. Thus, 



7_L<9 t 7^ , ''' + j±v x d x ir vv + j±v y d y ir 



nn 



' 7T 



TVV 
2 2 



rnv 

i 3 



7-L 



(A5) 



Note that in Ref. [29| this extra term was not present in 
Eq. (5.21a), but correctly added in Ref. The other 
relaxation equations, together with Iq, , , and l£ u 
and the vorticity tensor components remain formally un- 
changed. (This is so, since all nonvanishing Christoffel 
symbols are multiplied with u n = 0). 



APPENDIX B: (3+l)-DIMENSIONAL 
EXPANSION IN CARTESIAN COORDINATES 



This case is very similar to the two-dimensional case 
discussed in Sec. IIII CI The only difference is that now 
the velocity, spatial derivative, and all four-vector and 
tensor components in the z direction are non-zero. The 
velocity is = 7(1, v x , v y , v z ), where 7 = (1 — v 2 — v 2 — 

v z)^ 1 ^ 2 - Therefore, the new nonvanishing components of 
the charge four-current and energy momentum tensor, in 
addition to Eqs. ([M} - and Eqs. © - (fTU2"j) which 
formally remain the same with 7_l — > 7, are 



N z ee N°v z , 


(Bl) 


T Qz ee (e 


+ P) 7 2 « 2 + TT° Z , 


(B2) 


T zz ee (e 


+ P^vl + P + 7T" , 


(B3) 


T xz ee (e 


+ P)l 2 V x V z + TT XZ , 


(B4) 


T yz ee (e 


+ P)-f 2 V y V z + TT VZ . 


(B5) 


quantities 


are calculated similarly 


to the 



(2+l)-dimensional case, thus 



n = N»yjl-i%-i%-v*, 
e = (T oa - tt 00 ) - v x (T 0x - ir 0x ) 



(B6) 
(B7) 



Vy(T 0y 



While the velocity components in x and y directions re- 
main the same as in Sec. IIII CI the velocity component 
in z direction is 



T 0z _ n 



(B8) 



z yoo _ ^00 _i_ p ' 
The charge conservation equation is 

d t N° + d x {v x N°) + d y (v y N°) + d z (v z N°) = . (B9) 
The energy-momentum equations are 
8 t T ao + d x (v x T 00 ) + dy(v y T 00 ) + d z (v z T 0Q ) 

= -8 X ( Vx P - V X TT 00 + TT 0X ) - 8y (v y P - V y 7T°° + TT^) 

-d z (v z P - v z ir oa + ir 0z ) , (BIO) 
d t T 0x + d x (v x T 0x ) + dy(v y T 0x ) + d z (v z T 0x ) 
= -d x (P - v x ir 0x + ir xx ) - d y (-v y ir 0x + ir xy ) 

-d z (-v z tt 0x + n xz ) , (Bll) 
8 t T 0y + 8 x (v x T 0y ) + dy(v y T 0y ) + d z (v z T 0y ) 
= -d x (~v x n 0y + TT xy ) - d y (P - v y n 0y + 7T yv ) 

-0 Z (-v z TT ay + n yz ) , (B12) 
8 t T 0z + d x (v x T 0z ) + dy(v y T 0z ) + d z (v z T 0z ) 
= -d x (-v x n 0z + 7T XZ ) - d y {-v y TT°* + n yz ) 

~d z (P - v z tt 0z + n zz ) . (B13) 

The relaxation equations are formally similar to Eqs. 
(|119[) . (|120p . only the z-directed derivatives 7U z i9 z iI and 
r )v z d z ir^ 1 ' have to be added. Therefore, the new compo- 
nents of the shear tensor are 



J)z 



[yD(jv z ) + jv z D~f] + j 2 v z 



(B14) 



o" = -fl»(7«*)-7«*-D(7«*) + (l + r«j)g, ( B1 5) 

o xz = -\[d x {iv z ) + d z { 1Vx )] 

1 9 

- g h v xD(jv z ) + -iv z D(^v x )] + -i 2 v x v z - , (B16) 



) - V Z (T» Z - 7T U *) 



[d v (~fv z ) + d z (jv y )} 



- [yVyDfrvg) + -iVyD^v,)] + -i 2 v y v z - , (B17) 



where the expansion scalar is 6 = dt'y + d x (jv x ) + 
dy{l v y) + dz{iv z ), and the convective time derivative is 
D = u^d^ = •ydt + jv x d x + jv y d y + jv z d z . The form 
of the other components does not change in comparison 
with Eqs. (fl2T|) - (Tl^)) . 
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The term = (tt^u" + tt Xv u^)Du x leads to 

i°° = 2 1 [<K m D 1 -n Qx D( 1 v x )-iT 0y D{ 1 v y ) 

- ir 0z D( 1Vz )] , (B18) 
I\ x = 1 [(n 00 v x + ir^Dj - (ir 0x v x + ,«)%) 

- (77^^+^)13(7^) 

- (7r°^+ 7^)^(7^)] , (B19) 
i^ = 7 [(tt°% + tt ")^ - {n 0x v y + 7r xy )^(7^) 

- (tt°% + TT VV )D( 1 V y ) 

- (7r 0z % +7r^)£>(7^)] , (B20) 
i°* - 7 [(tt 0C V + 7r 0z )^ 7 - (tt ^ + 7r^)i?(7^) 

- (ir 0y v z + ir yz )D(jv y ) 

- (tt 0z v z + tt zz )D( 7 V z )] , (B21) 
If x = 2 1Vx [tt 0x Dj - tt xx D(jv x ) - n^D^Vy) 

- 7r xz D(jv z )} , (B22) 
If = 2-fVy [ir 0y D-/ - Tr xv D{ 1Vx ) - ix yy D(jv y ) 

- n yz D(jv z )] , (B23) 
if = 2 7 v z [n 0z D-f - Tr xz D{ 1Vx ) - n yz D{ 1Vy ) 

- ir zz D( 1Vz )} , (B24) 
if" = 7 [(7r°% + 7r% x )£> 7 - (tt*"^ + tt x V)£>(7^) 

- {ir xy v y + Tr yy v x )D^v y ) 

- (ir xz v y + ir y *v x )D(jv z )} , (B25) 
Jf* = 7 [(7r°^ z + tt '^)^ - (tt*^ + 7r a, *T; x )U(7« B ) 

- (tt^ + n vz v x )D(rfv v ) 

- (tt^ + tt^Mt^)] , (B26) 
if = 7 [(tt°^ 2 + 7r 0z v y )£»7 - (t^w, + n xz v y )D(~/v x ) 

- (ir yy v z +ir yz v v )D( 1 Vy) 

-- (7r^« z +7T zz % )£>( 7 ^)] . (B27) 



The terms Iq and 1%" are given by Eqs. (|2"Tjl and ([23ll . 
The new components which need to be computed com- 
pared to the previous case are Z^ 2 , Jf z , I yz , and if z . The 
non-vanishing components of the last term are 

7 3 00 = 2 (ir 0x u; x + ir 0y u; y + it 0z lu q z ) , (B28) 
I° x = n 00 uj x + n 0y u J x y + n 0z uj x z 

+ ir xx LU x + ir xy u y + ir xz uj z , (B29) 

4 y = 7r°V + n 0x u; y x + n 0z u; y z 

+ n xy uj x + n yy uj° y +7r yz LU z , (B30) 

i 3 ° Z = 7T °CJ 2 + TT° X U; Z X + TT° y LU Z y 

+ ir xz oj x + ir yz Lj° y + ir zz uj o z , (B31) 
I xx = 2 (n 0x uj x + n x y uj x y + tt xz u x z ) , (B32) 
If = 2(7r 0y u J v + 7r xy uj y x +Tr yz u J y z ) , (B33) 



i 3 = Z [TT U) q + 7T 



1.1 U IT (jj J- 7]-» [J 



7T Wj + 7T U! x + TT LO i) z 

^Qy , ,x 1 „yy , ,x , „yz , ,x 
7T "W o + 77""^ + 7T" U) z , 



(B34) 
(B35) 



/ 3 M = 7r 0x w z + 7r ax w z x + tt'^ 

+ TT 0z u x +ir yz uj x y +ir zz uj x z , (B36) 

If = Tr 0y u z +7r xy cj z x +TT yy cj z y 

+ n 0z Lu y + Tr xz u; y x + iT zz Lu y z , (B37) 

where the new vorticity tensor components arc 

1 



[dzl + dtiivz) + lv z D 1 - jD(jv z )] , (B38) 



u x z = 2 [dz(jv x ) -d x (jv z )] 

+ \ hv z D(jv x ) - yv x D(-yv z )] , 

U V Z = ^ [dz{lVy) ~ dy(rfV z )] 



(B39) 



(B40) 



such that uj° z = uj z — ~u 0z = uq z , uj x z = —uj z x — 
—oj xz = —oj xz and oj v z — —u z y — —uo yz = —uj yz - The 
other components are given in Eqs. (I140p - (|142p where 
one has to replace 7^ with 7. 



APPENDIX C: (3+l)-DIMENSIONAL 
EXPANSION IN (r,x,y,n) COORDINATES 

The metric of the space-time is the same as in App. 
lAl only the definition of the flow velocity changes. 
In this case the contravariant flow velocity is u M = 
r y{\,v x ,Vy,v 7] ), and the covariant flow velocity is it M = 
gfiuu" = 7(1, -v x , —v v , -t 2 v v ), where 7 = (1 — v% — v% - 
t 2 v 2 )~ x / 2 . The gradients are, — (d T , d x , d y , d v ) and 
d» = g^d v - (d T , -d x , -d v , -T- 2 d v ). 

Similarly as before the equations can be easily obtained 
from the ones found in Cartesian coordinates in App. [B] 
All laboratory frame quantities are formally the same, 
except for T zz -> T m = (e + P)j 2 v 2 + P/t 2 + tt"". The 
LRF quantities are 



n = N°^\-v 2 x ~v 2 y ~T 2 v 2 , 
e = (T m - 7T 00 ) - v x (T 0x - jr 0x ) 



(CI) 



- V y (T 0y -TT° V ) -T 2 V n (T° r > -7T°"). (C2) 

The velocity components v x and v y are given by Eqs. 
(|106[) . (|107p . the velocity component in 77-direction is 
given similarly as in Eq. (|B8[) . The charge conservation 
equation is given by 



d T N T + d x (v x N T ) + d y {v y N T ) + d v {v v N T ) 

= -V. 



(C3) 



The equation for energy-momentum conservation leads 
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to 



d T T TT + d x (v x T TT ) + d y (v y T TT ) + d n {v v T TT ) 
= -d x (v x P - v x n TT + ir TX ) - d y (v y P - v y ir TT + ir Ty ) 

-d n {v v P - v v tt tt + ^) - - (T TT + t 2 T^) , (C4) 
d T T TX + d x (v x T TX ) + d y (v y T TX ) + d v (v v T TX ) 

= -a x (p - v x tt tx + ir xx ) - d v (~v y 7T TX + n x y) 



-r,,*™ +***)- -17*, 



(C5) 



jy y u v 

d T T Ty + d x (v x T Ty ) + dy{v y T Ty ) + d v (v v T Ty ) 

= ~d x (-V X 1T TV + TT XV ) - d y {P - VyTT Ty + 7T yy ) 

-d v (-VrjTT^ + 7T OT ) - -T V , (C6) 

d T T™> + d x {v x T^) + d y {v y T^) + d v {v n T T7) ) 
= -d x (-v x 7r Tr > + n xr >) - d y (-vyir^ + n yr >) 



~9 V [ -J - V n TT Tri + TT 



Tin \ _ _jiTn 



(C7) 



The relaxation equations for the bulk viscous pressure 
and the shear stress tensor components ir xx , tt vv , ir xy are 
formally the same as in Cartesian coordinates, however, 
for the other components we obtain 



Dtt tt = 


-2r7^7r T " + I TT , 




(C8) 


Dn TX = 


~T~fv v ir xr > + r x , 




(C9) 


DTr Ty = 


-T~fv n ir m + r y , 




(CIO) 


Dn Tr < = 


-T^v v -K m - -ir Trl - 

T 


T 


(Cll) 


Dn xrl = 


7 'xri l v n tx | 
7T ' 7T + 

T T 


TXT] 

± 5 


(C12) 


Dn m = 


T T 


rvn 


(C13) 


Dir rm = 


-2l 7r 7)?) - 2^7r rr > 

T T 


+ r m . 


(C14) 



where denotes the right-hand side of Eq. (|T8|) , but in 
this case D — jd T + jv x d x + "yv y d y + 7i> ?) <9, 7 denotes the 
convective time derivative of scalars. 

The shear tensor components a xx , a vy , and a xy remain 
formally unchanged from Eqs. (fl~2"4"l) , (fl~25)) . (fl~26)) , while 
the ones which are different are calculated from Eq. ©, 



°™ = - T -^f^ + \[dAlv x )-d xl ] 

1 6 

- 2 bf D (T"*) + l v x D l\ + 7 2 ^3 : 

T^vtv,, 1 
V TV = -- L ^ UL + ^[dr{jV v )-d yl ) 

1 



TfJ 



T If 



1 9 

2 [7-0(7^) + lv v Dj] + 7 2 «r,g , 

-^(l + 2 7 ^r 2 )-i^( 7U ,) 
1 



..2„.2^ 



(C16) 



(C17) 



(C18) 



(C19) 



f 



1 6 

- [jv x D(~/v v ) + jv n D{jv x )} + -i 2 v x v n - , (C20) 

„3 



TVyVy _ 1 
T 2 



9 a ( 7 ^) + -2 d v(l v v) 



1 



- 2 [7«y-D(7^) + 7^)^(7%)] + 7 2 ^-fr, - , (C21) 

where the expansion scalar is 9 = 7 / V + (9*7 + 9 x (7fx) + 
9y('yvy) + dri('yv ri ). The Jo, /j"", ^2"; an d ^3"^ components 
remain formally the same. The new vorticity tensor com- 
ponents are 



= \ K7 + a T (T 2 7^)] 

+ i [rV^-T^rVr,)] , (C22) 

= J K(7^) - ^(T 2 7-y r) )] 

+ i [tV„.D(7^) - 7^i?(T 2 7^)] , (C23) 

= 5 K(7fy) -9 y (r 2 7^)] 

+ i [r 2 7 ^£>(7^y) ~ 7^£»(r 2 7^)] , (C24) 



a TT = -T 7 3 t; 2 + d r -f - 7Z37 + ( 7 2 - 1) - , (C15) where uj\ = u\, u x n = -u% and w» = -cj" y 
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